General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



NASA CR-1 44941 


SATELLITE LIFETIME ROUTINE USER’S MANUAL 

• by H. U. Everett and T. R. Myler 


(‘lASA-CF-l UUQUT) SV^LLITE LT7FTI °OT)TTNE N76-17172 

DATE'S HMIOAL (L^V Corn.) U) p HC 

*U.OO CSCT. 22? 


Unclas 

G 3/1 3 1U20S 



Prepared under Contract No. NAS1-12500 

VOUGHT SYSTEMS DIVISION 
LTV AEROSPACE CORPORATION 
Dallas, Texas 75222 

for 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


TABLE OF CONTENTS 


Page 


SUMMARY 1 

1.0 INTRODUCTION 3 

2.0 NOTATION 5 

3.0 PHYSICAL ENVIRONMENT * 7 

3.1 Earth Model 7 

3.1.1 Gravity 7 

3.1.2 Atmosphere 7 

3.2 Lunar Ephemeris 7 

3.3 Solar Ephemeris 9 

4.0 COORDINATE SYSTEMS 11 

4.1 Description 11 

4.2 Disturbing Acceleration Transformation 11 

5.0 PERTURBATION METHODS 13 

5.1 Lunar Disturbance 13 

5.2 Solar Disturbance 14 

5.2.1 Gravity 14 

5.2.2 Solar Radiation Pressure 15 

5.3 Earth Disturbance 15 

5.3.1 Gravity 15 

5.3.2 Atmosphere 17 

5.4 Total Derivatives of the Elements 18 

5.5 Numerical Integration Technique 18 

6.0 USER INSTRUCTIONS 19 

6.1 Introduction 19 

6.1.1 Limitations 19 

6.1.2 Computer Time Estimation 19 

6.1.3 Control Cards 19 

6.2 Input Data 20 

6.2.1 Rules for Data Card Preparation 20 

6.2.2 Table of Names 21 

6.2.3 Epoch Conditions 21 

6.2.4 F^rturbation Options 22 

6.2.5 Atmosphere Definition 22 


iii 


TABLE OF CONTENTS (Continued) 

Page 

6.2.6 Earth Model Definition 23 

6.2.7 Satellite Definition 23 

6.2.8 Termination 23 

6.2.9 Integration Controls 24 

6.2.10 Input and Output Control 24 

6.2.11 Error Exit Control 25 

6.3 Output Definitions 25 

6.4 Sample Problem 26 

APPENDIX A Atmosphere Model 33 

APPENDIX B Scientific Data Processing Routine Summary Documentation 39 

REFERENCES 41 


iv * 


SATELLITE LIFETIME ROUTINE 


USER'S MANUAL 


By H. U. Everett and T. R. Myler 
LTV Aerospace Corporation 


SUMMARY 


This report describes a FORTRAN coded computer program which determines secular 
variations in mean orbital elements of Earth satellites and the lifetime of the orbit. The dynamical 
model treats a point mass satellite subject to solar and lunar disturbing gravitational fields, second, 
third and fourth harmonics of the Earth's oblate potential. Earth's atmospheric drag and solar 
radiation pressure. Each of these disturbing functions may be selectively simulated. Data prepara- 
tion instructions, a sample problem and definitions of output quantities are included. 


1.0 INTRODUCTION 


This document presents a computing procedure for determining long period variations in the 
orbital elements of an Earth satellite. A time history ^ orbital elements yields the oroital lifetime 
when the perigee altitude decreases to near zero. Thus, one of the primary purposes of this computing 
procedure is to determine the lifetime of an Earth orbit 

Instantaneous derivatives of the orbital elements are calculated at equal time increments over 
one orbit period and averaged to obtain the secular derivatives which are integrated to yield mean 
orbital element time histories. The disturbing forces causing the orbital element changes are due to 
the sun, moon. Earth oblateness. Earth atmospheric drag and solar pressure. The secular derivatives 
of elements due to each of the disturbances, except solar pressure, are individually calculated and 
summed to obtain the total derivatives. 

A computer program similar to the procedure described herein was obtained from the NASA 
Goddard Space Flight Center and is described in Reference 1. The program of Reference 2 is an ex- 
tension and modification of the NASA version. Extensions were made to permit simulation of solar 
radiation pressure, third and fourth Earth oblateness harmonics and atmospheric drag. Modifications 
were to provide automatic stepsize and error control for numerical integration and to simplify pro- 
gram input and output. The current program is a modification of Reference 2 to remove computa- 
tional singularities in time rates of change of eccentricity and argument of perigee for circular orbits 
(zero eccentricity). 


Preceding page blank not ftlmhi 


2.0 NOTATION 


Symbols used in this report, excluding the appendices, are listed below with their definitions 
and units for program computations. 


A x ,A y ,A z 


C 

C D 

E 

e 

GM 

9 

H 

h 

I 

J 

K 

k 

4 

M 

P 

R 

R E 

r 

5 
t 
u 
V 
W 

w 

X, Y, Z 


atmospheric drag acceleration, R E da/“ 2 

Cartesian components of disturbing acceleration in the inertial equatorial frame, 
R e day" 2 

orbital semimajor axis, R E 

circumferential component of disturbing acceleration, R E day -2 
coefficient of drag 
eccentric anomaly, radians 
orbital eccentricity 

universal gravity constant, R E 3 day -2 
mean anomaly 

coefficient of third harmonic of Earth oblate potential 
e sin 03 

orbital inclination, radians 

coefficient of second harmonic of Earth oblate potential 
coefficient of fourth harmonic of Earth oblate potential 
e cos 03 

Solar mean longitude, radians 

mass of gravitating body, Earth mass M £ is the unit 

orbital semilatus rectum, R E , a ( 1 — e 2 ) 

c _2 
radial component of disturbing acceleration, R E day 

Earth equatorial radius, also a unit of length 

radius, R E 

aerodynamic reference area, m 2 
time since Jan. 0, 1961, days . 
argument of latitude, radians, v + 03 
inertial velocity, R E day 

component of disturbing acceleration directed normal to the orbital plane, 

R e day -2 

sine of satellite equatorial latitude 
Cartesian coordinates of position, R E 


IfcfiCJD 


'-Qvg 
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a 

P 

7 

e 

X 

P 

Q 

v 

a 

12 E 

cj 

Subscripts 

D 

E 

M 

0 

R 

S 

T 

X. Y,Z 
Superscripts 

F) 

( A ) 

( )' 


angle between geocentric radii to the satellite and to the sun, radians 
angle between geocentric adii to the satellite and to the tangent point on the Earth 
surface of a line from tf e satellite to the Earth-sun line, radians 
satellite illumination angle, radians 

obliquity of the ecliptic (angle between equatorial and ecliptic planes), radians 

satellite equatorial latitude, radians 

true longitude of the sun, i adians 

— 3 

atmospheric density, kg m 

2 —2 

oblate Earth gravity potential, R E day 
true anomaly, radians 

right ascension of the ascending node, radians 
inertial Earth rotation rate, 6.300383 radians day -1 
argument of perigee, radians 

quantity pertaining to atmospheric drag 
quantity pertaining to Earth 
quantity pertaining to moon 
quantity pertaining to Earth oblateness 
earth relative 

quantity pertaining to sun 
total 

component in coordinate direction 

derivative with respect to time 
unit vector 

components in ecliptic plane 



3.0 PHYSICAL ENVIRONMENT 


3.1 Earth Model 


3.1.1 Gravity. - Earth gravitational potential giving a disturbance to two-body motion is 


defined by terms through the fourth harmonic: 


<j> = GM, 


R, 


r 

+ K R 
30 r 2 


7 ? 


J (1-3 sin 2 X) - hr e (3— 5 sin 2 X) sinX 

3 5r 


(3.1) 


2 j 

E (3—30 sin X + 35 sin 


4 X) 


where 

Earth equatorial radius, R E = 6378.166 km 
Earth mass, M £ = 1 Earth mass ' 

Coefficient of 2nd harmonic, J = 0.00162345 
Coefficient of 3rd harmonic, H = -5.95 E— 6 
Coefficient of 4th harmonic, K = 7.95 E— 6 
Universal gravitation constant, GM E = 11467.849 R E 3 /day 2 
X is the latitude of the satellite position 

sin X =_?. 

r 


3.1.2 Atmosphere . — A static atmosphere model is used to compute the atmospheric density 
for altitudes below 120 kilometers. The calculations are based on an altitude- temperature profile 
that approximates the 1962 U. S. Standard Atmosphere Model. A dynamic atmosphere model is 
available for computing densities at altitudes above 120 kilometers. The dynamic model varies with 
time, location, and solar activity and is based on the 1969 NASA model presented in Reference (3). 
The computational algorithms used for both models are discussed in Appendix A. 


3.2 Lunar Ephemeris 

The lunar ephemeris is defined by the following expressions for mean elements presented 
in Reference (4): 

A m = 270.434358 + 13.1763965268d - 0.001 133T 2 + 0.0000019T 3 deg 
r M = 334. 329653 + 0. 1 1 1 4040803d - 0.01 0325T 2 - 0.0000 1 2T 3 deg 
= 259.183275 - 0.0529539222d + 0.002078T 2 + 0.000002T 3 deg 

where 

A m is mean longitude 

f M is mean longitude of perigee 




n M is longitude of mean ascending node 

d is ephemeris days from epoch of 1900 Jan. 0.5 Ephemeris Time (E. T.) 

T is Julian centuries of 36525 ephemeris days from epoch of 1900 Jan 0.5 E. T. 
M m = mean anomaly = A M — 

= true anomaly = M M + 2e M sin M M + ^e M 2 sin 2 M m + . . • 
co M = argument of perigee = r M — 

r M = raaius ■ P M 

1 +e M cosP M ! 


p M = semilatus rectum 

e M = eccentricity = 0.054900489 

Remaining constants for the moon or its ephemeris are 
semimajor axis, a M = 60.2681 R E 
inclination to ecliptic, sin l M = 0.089683448 
cos l M = 0.99597032 
mass, M m = M E /81.335 • 

Ecliptic Cartesian components of the moon position are 



where 


>M 


11 




12 


M 


13 


= cos(w M + t ; M )cosn M -sin(cj M + »> M ) sinft M cosl M 
= - sin (co M + j> m ) cos -cos (w M + v M ) sin cos l M 

= sin sinl M 


The equatorial Cartesian components are 


1 

X 

2 

1 


s 

ll 

y m 

li 

1 1 

V 

_ Z M_ 


_2m'_ 


(3.2) 
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where 



A 22 = cos e 
A 23 - -sin e 

A 31 =0 


Aj 2 3 sin 6 

A33 = COS € 

3.3 Solar Ephemeris 

Reference 5, page 98 gives the true longitude of the sun as a function of time (t) in days 
fram 1961 Jan. 0.0 as 

Xs = L s + 2e s sin gg 

where mean longitude is given by 

L s = -1.4062711 + 0.01 7202791 4t radians 

Solar mean anomaly is 

g s = -0.0496208 + 0.01 7201 9697t radians 

and solar orbital eccentricity (e s ) is 0.0167255. 

Instantaneous geocentric solar distance is given by 

1- e s 2 ^ 2 +3e s cosgg 

according to reference 6. Remaining constants pertaining to the sun or its ephemeris are 
semimajor axis, a s = 23454.708 R E 

argument of perigee, cu s = 4.923277 radians 

obliquity of ecliptic, e = 23.452294 -0.0130125T-0.00000164T 2 

+ 0.000000503T 3 degrees 

where T is Julian centuries of 36525 ephemeris 
days from the epoch of 1900 0.5 E. T. 

mass, Mg = 3.32951. 3M E 



* 


The equatorial coordinates of the sun are 
X s =* r $ cos Xg 

Y s = r s sin Xg cos e 
Z s « Tg sin Xg sin e 
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4.0 COORDINATE SYSTEMS 


4.1 Description 

Equations of motion are applied in an inertial equatorial frame as shown in Figure 1. The 
positive X axis is along the mean equinox of date, Z is collinear with the mean North polar axis, 
and Y is in the equatorial plane completing a right-hand orthogonal system. 


The lunar and solar mean elements are referenced to the ecliptic frame which is also illustrated 
in Figure 1. The X' axis is collinear with the X axis and the mean equincx of date, the Z' axis is 
normal to the ecliptic plane in the direction of the Earth's orbital angular momentum vector, and the 
Y' axis completes the right-hand orthogonal system. The angle between the equatorial and ecliptic 
planes is e as defined in Section 3.3. The coordinate transformation from ecliptic to equatorial 
components of a vector is given by 


X 

Y 

Z 




where [ A ] 


is defined in equation (3.2). 


4.2 Disturbing Acceleration Transformation 


Transformation of the inertial equatorial components of disturbing acceleration (A X ,A Y ,A Z ) 
to a satellite trajectory relative set is given by 


where 


R 


T 

m 

CM 

co 

CO 

1 


A X 

C 

= 

B 21 8 22 8 23 


a y 

W 


8 31 8 32 8 33 


A Z 


B 1 } = cos u cos ft — sin u sin ft cos I 
B 1 2 = — sin u cos ft — cos u sin ft cos I 
B 13 = sin ft sin I 

B 21 = cos u sin ft + sin u cos ft cos I 
B 22 = sin u sin ft + cos u cos ft cos I 
B 23 = - cos ft sin I 
B = sin u sin I 
B 32 = cos u sin I 
B 33 = cos I 


(4.2) 


Component R is along the radius, positive away from the origin; C is in the orbital planp and normal 
to R, positive in the direction of satellite motion; and W is normal to the orbital plane, positive in the 
direction of angular momentum. Accelerations R,C,W are based on the accelerations A X ,A Y ,A Z of 
individual disturbances which are defined in Section 5.0. 
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North Pole 


Z 



X, X' 

Mean equinox 

FIGURE 1. - LUNAR AND SOLAR ORBITAL GEOMETRY IN THE EQUATORIAL 
REFERENCE FRAME 
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5.0 PERTURBATION METHOOS 


The disturbing forces considered are those due to solar gravity, lunar gravity. Earth oblateness; 
Earth atmospheric drag and solar radiation pressure. The method used to simulate each of the dis- 
turbances is a technique of averaging the orbital element derivatives over an orbit period at equal inter- 
vals of mean anomaly. This method accounts for the inertial position of the moon and sun according 
to the date and time. Additionally, this method permits simulation of solar radiation pressure. 

The classical set of orbital elements normally used for analyses include: a (semimajor axis), 
e (eccentricity), go (argument of perigee), £2 (right ascension of ascending node) and I (inclination). 
However, because the derivatives of e and go are undefined at zero eccentricity, parameters h and k 
as described in Reference (7) are .lumerically integrated in the computational procedure. These 
parameters are defined in terms of e and go as foilows: 

h = e sin go 
k = e cos go 

Derivatives of h and k are taken from the above mentioned technical paper and are not a function of 
e or go. Since h and k are integrated, e and go are obtained from the following relationships 

e = V h 2 + k 2 
go = tan -1 

Derivatives of these quantities which are required for program output are 

. hh +kk 

e = 

e 

kh — hk 



5.1 Lunar Disturbance 

Inertial accelerations on the satellite due to the lunar gravity are as follows: 


A x = — GM 


U M 3 r M / 


Y.Z 


(5.1) 


where 


V = < X - X M> J + 


I 
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These inertial accelerations are transformed to the satellite trajectory relative set using equation (4.2). 
Then the instantaneous derivatives for the orbital elements are obtained from R, C, W as shown 
below. Derivatives for a, ft and I are taken from Reference 8 and derivatives for h and k are taken • 
from Reference 7. 



Averaging the instantaneous derivatives over the orbit period is done in equal increments of mean 
anomaly; and therefore equal time intervals. The secular rate, averaged with respect to mean 
anomaly, is 

^ SEC = J_ f 2?r ^2. dg ft — ► I, a, h,k 

2jt J dt 
0 

where g is mean anomaly. A transformation from mean to eccentric anomaly is made to avoid 
repeated solution of Kepler’s equation for mean anomaly. 

Since g = E - e sin E 

then dg = (1 — e cos E) dE dE 

a 

therefore . 2ir 

ft$gc = J_ T r dft dE ft — I, a, h, k (5.3) 

2rra J dt 
0 


5.2 Solar Disturbance 

5.2.1 Gravity. — Inertia; accelerations of the satellite due to solar gravity are as follows; 


A x GM S 


£?: 


T) 


X-Xc 


V 


+ _^S_ 1 X — ► Y, 

rc 3 J 
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where 


Ag 2 = (X-X s ) 2 +( Y-Y s ) 2 + (Z-Zg ) 2 

SFACT is solar radiation multiplier defined in Section 5.2.2. 


These inertial accelerations are transformed to the satellite trajectory relative set using equation (4.2). 
Instantaneous derivatives of the orbital elements are obtained using equations (5.2) and the secular 
derivatives are obtained using equation (5.3). 


5.2.2 Solar Radiation Pressure. — Solar radiation pressure is simulated by reducing the GMg 
when calculating the acceleration on the satellite due to the sun. During satellite, illumination 

SFACT = 1 — fiM. 

GM e 


Sk , 2 

GM= "w ^. 9 ° 


S = satellite reference area towards sun in m 2 

k = solar flux constant at 1 AU, 

1.03C34x 10-6 lbf/m 2 

W = satellite mass in lb. 

3g = semi-major axis of solar orbit in R E 

9 0 = acceleration of gravity at Earth's surface in R E /day 2 

When the satellite is in the Earth's shadow, SFACT is defined to be unity. Figure 2 shows the geome- 
try used to determine when the satellite is illuminated. Application of plane geometry laws gives 


7 = -l-a + 0 
2 

where A A 

cos a = r r s 


cos P = — £■ 
r 


When the satellite illumination angle 7 is negative, the sun is eclipsed by the Earth. Conversely, the 
satellite is illuminated when 7 is positive. 

5.3 Earth Disturbance 

5.3.1 Gravity. — Disturbing acceleration components in the inertial equatorial frame are 
obtained by differentiating equation (3.1) 

d<t> 

A* = — X Y, Z 

X dX 


15 


Sun 



FIGURE 2. - GEOMETRY TO DETERMINE SATELLITE ILLUMINATION 
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or 


A x =-X [F(Z,r)] 

A y = — Y [F(Z,r)J 

-gm_ e ^2 r-j w(3 _ 5w2) + h_ Re (30w 2 _ 35w 4 _ 3) 
1 W (15-70w 2 + 63w 4 )] 


A Z = 


where 


w = 


[f (Z, r) j - |\j ( 1 _ 5 w 2, + h «e (3 - 7w 2 )w 

+ (3 - 42w 2 + 63w 4 ) J 

The above inertial accelerations are transformed to the satellite trajectory relative set using equation 
(4.2). The instantaneous derivatives of the orbital elements are then obtained from equation (5.2). 
The secula. derivatives are obtained using equation (5.3). 

5.3.2 Atmosphere. - The satellite trajectory relative acceleration components due to atmos- 
pheric drag are obtained from 


R = — 


-A n e sin v 


^R ~\J 1 + e 2 + 2e cos v 
M ( 1 + e cos i>) 


-A n r v( 1+e 


2e cos v 


— rSl E cos * 


-Ar 


W= — 2. n £ r cos u sin I 


where 


a d " 


V = 


1 pV R C Q S 
~2 


m 


-FHFT) 


V R = V - n E r cos I 
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The instantaneous derivatives of the orbital elements are then obtained from equation (5.2). The 
secular derivatives are obtained using equation (5.3). Atmospheric effects are neglected if the satellite 
is above 1400 kilometers. 


5.4 Total Derivatives of the Elements 

Total derivatives for the orbit elements are obtained by simply summing contributions from 
the various sources of disturbance. 

*T = a D + a S + a O + a M 

*T = *D + *S + 'o + *M 
• • • • • 

+ J2q + n M 

• • • • • 
hj = h D + h s +h Q + h M 

• • • • • 

In most instances, the derivatives a s , a Q and a^ should vanish except for numerical inaccuracies and, 
therefore, are not usually used in the equation for a-p. However, the a s and a M may be selectively 
used if desired ty the user. 


5.5 Numerical Integration Technique 

A Runge-Kutta technique which includes error control and automatic interval sizing is used 
to numerically solve the deferential equations for the orbital element secular rates. Time is the 
independent variable. 

Development of the integration procedure occurred at the NASA Lewis Research Center and 
is described in Reference 9, Appendix D. 
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6 0 USER INSTRUCTIONS 


6.1 Introduction 

6.1.1 Limitations . - Neither □ ground track nor a time history of satellite position in any 
inertial frame of reference may be obtained since satellite position in orbit is not integrated. 

Mean orbital elements are used Ly the procedure; thus, short period variations in the elements 
having a frequency on the order of the satellite orbital period are not simulated. 

Since the perturbation method uses the averaging technique, care should be taken that a 
sufficient number of samples are specified to yield an accurate average. Required number of sampl- 
ing points will vary for different prob!ems. One test for accuracy is to jompare results from the 
same problem but with different numbers of samples per orbit. Probably a more satisfactory check 
is simply to verify the derivatives of semi-major axis that are smali enough tc be neglected— e.g., 
those due to solar and lunar disturbances for a near Earth orbit and that due to the oblate dis- 
turbance. 

Inclination must not be zero during the calculations since inclination is in the denominator 
of the derivative of ascending node. The minimum value of inclination has not been determined but 
an input value of 0.01 degrees has been used with success. An increased number of integration steps 
results from near-zero inclination due to the large derivative of ascending node. 

6.1.2 Computer Time Estimation . — Average computation time on the CDC 6600 is about 
six integration intervals per CP second. This time is applicable when simulating all four disturbances 
averaging at 30 points per orbit and includes the re-entry phase. For cases which do not include 
the re-entry phase, the computation time is about twelve integration intervals per CP second. 

6.1.3 Control Cards. — The control cards and deck setup required to execute two problems on 
the CDC 6600 computer are shown below. Extension to three or more problems is straight-forward. 

JOB, CM 115000, T200. 

PROJECT 

ATTACH (GO, S7051B, ID=EVERETT) 

GO. 

7/8/9 multi-punched in column I 

Table of Names 

$D=1 

Data Cards of Problem Set 1 

$D=1 

Data Cards of Problem Set 2 

$D=1 

6/7/8/9 multi punched in column I 



6.2 Input Data 


Program input is achieved by using the NASA input subroutine, reference 10, which utilizes 
arithmetic input statements. This subroutine provides flexibility and ease in inputting data. Rules 
for preparaing the input cards are shown below. 


6.2.1 Rules for Data Card Preparation . — 

(1) Input data have preass : gned names and storage locations. Data arc- input by a 
statement of the type: 

ISTART = -1, START (2) = 0, , 127.1/15, - 30E - 2 $$ SCOUT 


This card will store — 1 in the location identified as ISTART, zero in the second 
location of the START array, will not disturb the third location of the START 
array, will store the quotient of the indicated division in the fourth location of 
the START array and —0.30 in the fifth location of the START array. The $$ 
causes the card to print on the output listing. The word SCOUT appears as comment 
only. 

(2) The data card format is flexible. Blanks are ignored, except in the alphanumeric 
field. All 80 columns may be used. Decimal points are optional. 


(3) If the $$ is omitted, the data is stored but the card is not printed with the output 
Comments may be placed to the right of the $$, however, avoid any characters 
adjacent to the $$ since these may result in undesirable printer line spacing. 

(4) A comma after the last value on a card is optional if the next card begins with a 
variable name. 


(5) The following arithmetic operations are permitted: 


Addition 

Subtraction 

Multiplication 

Division 


Use the + character 
Use the — character 
Use the * character 
Use the / character 


Parentheses to indicate order of arithmetic operations are not permitted. The 
order of operations is from left to right 

(6) As indicated in the example data card, each variable may be regarded as an array. 

(7) Where built-in (BN) values are indicated in the input definitions, the parameter 
need not be input unless different values are required. All parameters are 
initialized zero unless a BN value is indicated. 

(8) Preceding and following each problem set must be a $D=1 card. A single $D=1 
card must separate problem sets. 
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(9) Alphanumeric or title information is input as follows: 

LINE 1 =(A1 4) F-1 TRAJECTORY$$ 

Fourteen is the number of columns containing the title, excluding the $$. 

(10) Units are kilometers, days, degiees, kilograms unless otherwise noted. 

(11) Preceding the $D=1 card of the first problem set must be the table of names. 

6.2.2 Table of Names. — The table of names assigns an input variable name to a single or 
several storage locations, which are assigned to variables internal to the computer program. All 
variable names appearing on the input data cards must be assigned a storage location in the table of 
names. The table of names is shown in Section 6.4 in the sample problem. 


6.2.3 Epoch Conditions. — 


JDATE (1) 

Julian Date of the epoch, ends in 0.5 

JDATE (2) 

fractional day of the epoch 

ISTART 

= —1, mean orbital elements at the epoch JDATE 
determined from the START array. 

=1, is the same as ISTART =-1 except START (4) 
is satellite right ascension 

=0, Mean orbital elements at the epoch are input 
directly as E, 1, ARGP, NODE and A or N 

START (1) 

radius to satellite, feet 

START (2) 

inertial velocity of satellite, ft/sec 

START (3) 

geocentric latitude of satellite position 

START (4) 

Greenwich East longitude of satellite position 

START (5) 

heading of satellite inertial velocity 

START (6) 

flight path angle of satellite inertial velocity 

A 

satellite orbital semimajor axis 

E 

satellite orbital eccentricity 

1 

satellite orbital inclination, must not be zero 

ARGP 

satellite argument of perigee 

NODE 

right ascension of ascending node of the 
satellite orbit on the equatorial plane 

N 

satellite mean motion, rev/day, used if 
A=0 and ISTART = 0. 
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6.2.4 Perturbation Options. — 


MOON 


SUN 


OBLATE 


DRAG 


NSUN 

SPRESS 


=0, lunar disturbance is ignored 

=1, uses a procedure which averages the lunar disturbance 
over the satellite orbit at equal intervals m satellite 
mean anomaly. ( 1 BN) 

=2, same as MOON = 1 except lunar disturbance effect 
on semi-major axis is included. 

=0, solar disturbance is ignored. 

=1 , uses a procedure which averages the solar disturbance 
over the satellite orbit at equal intervals in satellite 
mean anomaly. (1BN) 

=2, same as SUN=1 except solar disturbance affect on semi- 
major axis is included. 

=0, earth oblateness disturbance is ignored. 

=1, second, third and fourth harmonics of the earth's 
oblate potential are simulated by averaging the 
disturbance over the satellite orbit at equal intervals 
in satellite mean anomaly. (1BN) 

=0, earth atmospheric drag is ignored. 

=1, uses a procedure which averages the drag disturbance 
over the satellite orbit at equal intervals in satellite 
mean anomaly. (1BN) 

integer number of equal increments in mean anomaly for 
averaging disturbances over the satellite orbit to 
obtain secular rates. 1(KNSUN<360 

satellite reference area used in solar pressure calculation, 
square meters. SPRESS is the satellite area per- 
pendicular to the satellite - sun line. (1BN) 


6.2.5 Atmosphere Definition. — 

LUIGI specif ies atmosphere definition. Cannot be changed in 

subsequent cases. 

=0, atmospheric quantities are calculated based upon 
1962 Standard Atmosphere. 

= 1 or 2, atmospheric quantities are calculated based 
on an exospheric temperature profile based on 
local time, season and time within the eleven year 
solar cycle. Table of FBAR must be supplied. 

= -1, oi -2, atmospheric quantities calculated from an 
input constant value of exospneric temperature, 
TEXCON. 
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=±2, a table of density as a function of exospheric 
temperature and altitude is printed. 

TEXCON constant exospheric temperature, used when LUIGI a —1 

or —2, must be between 650 and 2100 degrees Kelvin. 

FBAR table of solar flux (F. j0 7 ), 100 values maximum, units 

are 10“ 2 2 watts/in 2 /cycle/sec. 

TCYCLE table of time for FBAR, units are years A. D. 

i 

6.2.6 Earth Model Definition. — 

REQUAT earth equatorial radius (6378.166 BN) 

EESQRO eccentricity squared of spheroidal earth model 

(0.0066934217 BN) 

OBLATJ coefficient of the second gravitational harmonic 

(0.00162345 BN) ! 

OB LATH coefficient of the third gravitational harmonic 

(-5.75E -6 BN) 

i 

OBLATK coefficient of the fourth gravitational harmonic 

(7.95E-6BN) ! . 

GM gravitational constant, ft 3 /sec 2 

(1.4076576E16 BN) 

6.2.7 Satellite Definition. - 

CD aerodynamic drag coefficient (2.5 BN) 

SREF aerodynamic reference area, square meters (1 BN) 

MASS satellite mass (1 BN) 

6.2.8 Termination. — 

TSTOP time since JDATE at which case will terminate. 

ITERM selective parameter for terminating the case when 

FINALVALUE is reached. Case will terminate 

on earliest of TSTOP, FINALVALUE and STEPMX. ! 

Dependent Variable 
=1 semi-major axis 

=2 eccentricity 

=3 inclination 

=4 argument of perigee 

=5 right ascension of ascending node 

=6 perigee altitude 

=7 time derivative of semi-major axis due to 

atmospheric drag, km/day. 
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FINALVALUE value of dependent variable for case termination 
when ITERM is non-zero. 

SLOPE = -1, termination occurs when ITERM parameter is 

decreasing. 

=0, termination occurs first time ITERM parameter is 
attained. 

= 1, termination occurs when ITERM parameter is 
increasing. 

XTOL tolerance on the independent variable, time, that * 

will cause case termination on a dependent variable 
when ITERM is non-zero. (0.0001 BN) 

6.2.9 Integration Controls. — 

BACKUP =1, causes the integration to proceed backwards in 

time. All inputs remain positive. 

EREF reference'value of normalized truncation error. 

(IE -4 BN) 

ERRFAC factor by which the normalized truncation error may 

exceed EREF before rejection of the interval occurs. 

<5 bn) ; 

DSTART Estimated initial integration interval. The integration 

interval is subsequently varied by the program to 
control truncation error. (1 BN) 

6.2. 1 0 Input arid Output Control. — 

NPROB problem number printed with page headings. It is 

incremented by the program when successive problems are 
run (1 BN) 

LINE 1 =(A67) Line one title in columns 12 through 78. 

(blanks BN) 

LINE2=(A67) Line two title in columns 12 through 78. 

(blanks BN) 

SAVE =1, causes all input data to be saved for possible 

subsequent problems. If SAVE=0 all input 
locations are cleared and restored to their 
built-in values between problems. 

DPRINT print interval for time. If DPR INT=0, output will 

occur every nth successful interval where STEPS=n. 

Abnormal termination may occur with large values 
of DPRINT. 
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number of integration steps between printouts 
when DPR INT=0. (1 BN) 

CHKOUT =1, prints out h, k, h, k parameters (integration 

variables) as a function of time, on page C. 

• • 

=2, in addition to h, k, h, k parameters, printout 
moon ephemeris (right ascension of ascending 
node and ueclination) at each time step. 

6.2. 1 1 Error Exit Control. - 

DMIN minimum continuous integration interval size. (0.5 BN) 

STEPMX maximum number of successful and unsuccessful integration 

intervals which may occur during one problem. (100 BN) 

ESTART =0, causes the program to dump COMMON data and 

• exit the computer when an error is encountered 

=1, causes the program to du np COMMON data and proceed 
to the next problem when an error is detected. 

= —1, causes the program proceed directly to the next 
problem when in er. or is detected. 


Input Options and Data 
SFACT 


Page A 
TIME 
ARG PER 

ASC NODE 

INCL 
ECCEN 
SEMIMJR 
PER ALT 


6.3 Output Definitions 
Self Explanatory 

1. — GM/GM S where GM = solar gravitational 

constant acting on satellite in 

sunlight. See Section 5.2.2 

GM S = unive r sal gravitational constant 

time from the epoch, days. 

satellite argument of perigee, deg. (zero is printed whan 
eccentricity is zero) 

right ascension of the satellite ascending node in the 
equatorial plane, deg. 

inclination of the satellite orbit to the equatorial plane, deg. 

eccentricity of the satellite orbit 

semimajor axis of the satellite orbit, km. 

perigee altitude of the satellite orbit above a sphericai 
Earth, km. 
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STEPS GOOD 
STEPS BAD 


indicates which orbital element hat tne largest integration 
error 


NERR 

Element 


2 

semimajor axis 


3 

h parameter (e sin 

w) 

4 

inclination 


5 

k parameter (e cos 

w) 

6 

node 



count of successful integration intervals 
count of unsuccessful integration intervals 


SHADOW POINTS count of NSUN points in the current orbit which are 
eclipsed by the Earth 

DERIVATIVES OF SEMIMAJOR AXIS IN KM/DAY 
LUNAR lunar contribution to a, km/day 

SOLAR solar contribution to a, km/day 

OBLATE oblate Earth contribution to a, km/day 

DRAG atmospheric drag contribution to a, km/day 


Page B Page B output is self explanatory 


Page C Page C output is self explanatory. PageC 

applies to parameters h and k and is printed only 
if CHKOUT^O. 




6.4 Sample Problem 

Output of a typical problem is presented in this section. The problem is for an initial orbit 
of 250 km perigee, 1200 km apogee and an epoch of 20 November 1975. All orbit perturbations are 
considered and a solar flux time history is predicted for the atmospheric drag calculations. 

The first page is the table of names, which precedes the data of the first problem set The 
second page is a list of the data cards of the first problem set. These first two pages are printed since 
each card has the $S characters. Following these pages is a page showing input options and data, 
followed by Pages A, B, and C, as defined in Section 6.3. 
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APPENDIX A 


ATMOSPHERE MODEL 


Static Atmosphere Model 

A static atmosphere model is used to calculate the density for altitudes below 120 kilometers. 
The atmosphere is defined by an assumed relationship between the temperature and the geopotential 
altitude and by a sea level pressure. Geopotential altitude is the altitude above a constant gravity 
earth which gives the same potential energy as altitude above an earth with inverse square gravity. 
Geopotential altitude is 



where R Q is the mean earth radius and h is the geometric altitude. A continuous function with linear 
segments is used as the relationship between the temperature and geopotential altitude. For each 
segment 

T = T r + L n< H - H n > 

where L n is the slope of the linear segment and T n is the temperature at altitude H n . 

If the perfect gas equation of state, 

P-5l£ 

Mo 

and the hydrostatic equation. 


dP = -pg Q dH 


and the linear temperature-altitude relation are combined and integrated, the pressure is expressed as 


P = P n 
P = P n 



The density is 


P = 


™o 

RT 


L n =£0 

L n =0 


and the speed of sound is 
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The following constants are used to approximate the 1962 U. S. Standard Atmosphere. 


Sea Level Pressure 

101325 newtons/m 

Units constant, g Q 

9.80665 m/sec 2 

Sea Level Molecular Weight, M Q 

28.9644 

Universal gas constant, R 

8314.32 joules/ kg° 

Mean earth Radius, R Q 

6,356,766 m 

Specific heat ratio, 7 

1.4 

Geopotential 

Temperature 

Altitude, km 

°K 

0 

288.15 

11 

216.65 

20 

216.65 

32 

228.65 

47 

270.65 

52 

270.65 

61 

252.65 

79 

180.65 

88.743 

180.65 

98.451 

210.65 

108.129 

260.65 

117.776 

360.65 

146.541 

960.65 

156.071 

1110.65 

165.571 

1210.65 

184.485 

1350.65 

221.967 

1550.65 

286.476 

1830.65 

376.312 

2160.65 

463.526 

2420.65 

548.230 

2590.65 

630.530 

2700.65 
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Dynamic Atmosphere Model 

Because of the dynamic nature of the upper atmosphere and the importance of the effect of * 
aerodynamic drag on orbital motion, a dynamic atmosphere model is used for altitudes above 120 
kilometers. The model is a slightly modified version of the 1969 NASA model. Reference (3), which 
is based on the model of Jacchia, Reference (11). The modifications are as recommended in Refer- 
ence (12). The model requires an input table of 10.7 cm mean solar flux and a table of geomagnetic 
index. The computational algorithm, which is almost identical to that of Appendix A of Reference 
(3), is presented below without discussion. 


A. Exospheric Temperature Computation 

1. Angle between atmospheric bulge and computation point 

t = H - 45° + 1 2° sin (H + 45°) (± 180°) 

where 

H is the hour angle of the computation point 

2. Mean solar activity correction. 

T 1 = 362 + 3.60 (?) 

where F is the input value of mean solar flux. 

3. The daily solar activity correction is neglected. 


4. Semi-annual correction 


T 3 = T 2+ f F 


where 
f = 


\ 0. 37 + 0. 


14 sin 


27T 


P— 151 \ ) 
365 J ) 


)]j *[“ ia 


D is day number. 

5. Diurnal Correction 


T 4 = T 3 


1 + 0.28 sin 2 - 4 5 9 

1 + A cos 2 5 (l~\ 


\ 2 / _ 
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where 
A = 0.28 


Fcos^ 1 ® W — sin 2,5 0~ 1 
|_1 + 0.28 sin2-5 g j 


W=— (X-6) 
2 


0 


1 (X + 5) 

2 


X = latitude of computation point 
5 = declinat : on of the sun 


6. Geomagnetic activity correction 

T 5 = T 4 + La p + 100 £l - exp (-0.08a p )J 
where 

L = 1 + 2.85 (X - 30°) for |x|>30° 
L = 1 for | X | < 30° 

a p = input value of geomagnetic index. 

B. Temperature at given geometric altitude 

T 6 = T 5 - [t 5 - 355] [exp l-SAH)] 
where 

AH , (Z— 120) (6476.7 7) 

6356.77 + Z 


Z = geometric altitude, km. 

S= 1.5 x i0- 4 + 0.029 exp (— X 2 /2) 

T 5 - 800 

X= 750+ 1.722 x 10 ~ 4 (T 5 -800K 
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C. Number Density Computations 

1. Thermal diffusion factor for hydrogen 
T 0 = -10.48947 + 2.844291 x10" 2 [ 


- 3.620958 x 10“ 5 [j 5 J 2 + 2.341 193 x 10“ 8 T g 3 

- 7.577509 x 10“ 1 2 £t 5 ] 4 + 9.753963 x 10” 1 6 T 5 " 5 


2. Hydrogen number density at 500 km altitude 


= antilog 10 J^73. 13 - 39.4 log 1Q T & +5.5{log 10 T 5 ) 2 J 


3. Hydrogen number density for altitudes greater than 500 km 

, , , f K1 + T n + 1.008Q) r 1 

N(H) = N(H) 500 [bJ d [exp(-1.008SAHQ)J 


where 


1 - P 

1 — p exp (-SAH) 


T 5 - 355 


1.13619 


4. Helium number density 


7 T “I (0.63 + 4.002Q) T "I 

N(HE) = (3.4 x 1 0 7 ) [bJ I exp(-4.002S AHQ)J 


37 


1 I I I i I l 

APPENDIX A I 

I 

i 

I 

5. Number density for molecular nitrogen and molecular and atomic oxygen. 

N(l) = [n(I) 120 ] [b] [ 1 + QM(I) ] jexp[-2AHQM(l)]j 

I n 2 , o 2 , 0 

where 

| 

N( N 2 ).2o =40x 1011 cm' 3 i 

N(° 2 )i 2 0 = 7 - 5x 10 10 cm' 3 

N(0>120 = 7.6 x 10 10 cm’ 3 

M(N 2 ) = 28.0134 

M(0 2 ) = 31.9988 

M(O) = 15.9990 

l 

D. Mass Density 

p = N(H) W(H) + N(HE) W(HE) + N(N 2 ) W(N 2 ) | 

+ Ni0 2 ) W'0 2 ) + N(O) W(0) gm/cm 3 

where I 

W(H) = 1.6731 x 10" 24 gm/mole 
W(HE) = 6.6435 x 10" 24 gm/mole 
W(N 2 ) = 4.6496 x 10" 23 gm/mole 
W(0 2 ) = 5.3104 x 10" 23 gm/mole 
W(0) = 2.6552 x 10" 23 gm/mole 
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SCIENTIFIC DATA PROCESSING ROUTINE 
SUMMARY DOCUMENTATION 


IDENTIFICATION 


Title Satellite Lifetime Routine 

Routine No. 7051 Date Filed 1965 

Responsible Engineer H. U. Everett Unit 

Date Completed November 1975 Source 

Language: 


Security Class. U 


2-53012 Telephone Ext. 7694 


X FORTRAN Other 


Key Words orbit lifetime, orbit element derivatives, atmosphere 
RESOURCE REQUIREMENTS 

Typical CPU + I/O Time 20 sec Machine(s) CPC 660 0 


Core 115k oct al Tape 0 | 1 Plot | | Graphics 

DESCRIPTION (Include: Purpose of routine , i nput , output , and fu nctional 
description ) 

Purpos? — The primary purpose of the routine is to calculate time histories of the 
orbital elements and the orbit lifetime. 

input Description — Input data cards define the initial orbital elements or injection conditions, the 
ballistic coefficient of the satellite, the perturbing forces to be considered, and parameters 
of the dynamic atmosphere model. 

Output Description - Output consists of columnar time histories of the orbital elements and their 
time rates of change. 

Functional Description — The routine numerically integrates the secular (averaged) time rate of 
clianqe of the orbital elements to obtain their time history. The math model includes 
perturbing forces due to the sun, moon, oblate earth giavity, and a dynamic atmosphere. 
Atmospheric properties are influenced by solar radiation seasonal variations in atmospher- 
ic properties and the "bulge" in the atmosphere that is oriented with respect to the 
earth sun line. 

DOCUMENTATION 


No. 

Source Cards 2700 
Other 


VSD Report No. 2 53010/5R 23029, "Satellite Lifetime Routine User's Manual," 
dated 1 December 1975. 
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